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ABSTRACT 

Perturbation theory is used to calculate the frequency shift of acoustic 
modes of an homogeneous turbulent fluid and the frequency shifts of solar 
modes due to turbulent convection. For sound waves in a random flow, the 
fractional frequency shift is +(11/30)M 2 in the long-wavelength limit, and the 
shift at short wavelengths is — (2/3) M 2 , where M is the average Mach number 
of the flow. In our model of the solar convection, the low-degree f-mode shift is 
positive and unobservably small, whereas the fractional frequency shift at high 
degrees (I 3> 500) is — M 2 , where M is the Mach number of convection near the 
solar surface. 

Subject headings: turbulence — Sumoscillations 
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1. INTRODUCTION 

Frequencies of low-order (n ~ 1) high-degree (/ ~ 1000) solar modes are measured with 
an accuracy of about 0.2% (Libbrecht, Woodard & Kaufmann 1990, Bachmann, Duval, 
Harvey, & Hill 1995). It is found that the observed frequencies of the fundamental mode 
are smaller than the theoretically expected value of ^fgk by ~ 1%. One may expect that 
turbulent convective flows with characteristic velocities v ~ lkm/ s, in the region with a 
characteristic sound speed c ~ 10km/ s, will cause a fractional frequency shift of the order 
v 2 /c 2 , ~ 1% shift. Magnetic fields seem to be less important than turbulent motions, 
because they reach equipartition strengths only in sunspots. Also, frequency shifts due to 
magnetic fields should be positive. We will ignore magnetic fields in what follows. 

The negative sign of the shifts was explained by the following simple argument 
(Murawski & Roberts 1993; Rosenthal 1996). Consider a sound wave propagation through 
a layer of thickness 2L. Let the flow velocity along the wavevector be +V in the first 
sub-layer of thickness L, and —V in the second sub-layer. Then the total propagation time 
t = L/(c + V) + L/(c - V), and the effective sound speed is c* = 2L/t = (c 2 - 2V 2 ) 1/2 . 
Therefore the fractional frequency shift, s = Scu/cu, is 

S = ~\ M ^ (!) 

where M =< V 2 > /c 2 , and 1/3 comes from the angle average. 

The above kinematic explanation is incomplete. The sound speed in a random 
homogeneous isotropic flow depends on the wavelength. We show in §2 that the effective 
sound speed is increased in the long wavelength limit and decreased in the short wavelength 
limit. We discuss the solar frequency shifts in §3. Our perturbation theory (PT) calculations 
are in principle straightforward, but algebraically tedious; appendices explain how we 
performed the calculations and tested the results. 
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2. SOUND SPEED IN A HOMOGENEOUS ISOTROPIC RANDOM FLOW 

We will assume that the turbulence Mach number is so small that the turbulent flow 
is approximately stationary (even at time scales of the long wavelength acoustic modes). 
Only the turbulent velocity field, V, can give frequency shifts that are quadratic in M. 
Turbulent pressure, density and entropy fluctuations are quadratic in M, and give rise to 
0(M 4 ) frequency shifts that can be neglected. To within this accuracy, V • V = in the 
absence of stratification. 

In the leading order, one can find the shift caused by a single Fourier harmonic of the 
flow, and then average the result over all harmonics. Since velocity gradient of a Fourier 
harmonic of the incompressible flow is perpendicular to the velocity, it is sufficient to solve 
the problem for a shear flow 



The sound wave frequency in an isotropic random flow is then obtained by taking an 
appropriate average. We will use the fact that V is a pure Fourier harmonic later in our 
calculation, but at this stage V may be thought of as an arbitrary shear flow with a zero 
mean. 



We first obtain a convenient representation for the linear sound wave equation in the 
presence of a shear flow. Let the velocity perturbations (u, v, w) and dimensionless density 
perturbation (p) associated with the sound wave depend on x in an unspecified way, and be 
proportional to exp(— iut + ik y y + ik z z). Then, denoting the x-derivative by a prime, and 
using the standard notation fl(x) = uj — k y V(x), we have 




(2) 



2.1. The Basic Equation 



— iVtu = —c 2 p', 



(3) 
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- iVtv + V'u = -ic 2 k y p, (4) 

— iVLw = —ic 2 k z p, (5) 

— iVLp + u' + ikyV + ik z w = 0. (6) 

Equations (2)- (5) give the velocity perturbations in terms of the density perturbation. 
Using these expressions in Eq. (6), we get 

(^-kl) P + p"- 2 -^p' = 0, (7) 

where k\ = k 2 . + k 2 z . Up to the second order in V, this equation can be written as 

-p" + Q lP + Q 2P =( l 4-k 2 ± )p, (8) 



where Q\ is the linear perturbation operator 

2k 



Qi = ^(k 2 V-V'd), (9) 

and Q2 is the quadratic perturbation operator 

Q 2 = -^(k 2 V 2 + 2VV'd). (10) 

In the above expressions we have replaced uj by the unperturbed value ck, knowing that the 
frequency shift due to the linear perturbation Q\ vanishes in the first order PT. 



2.2. Perturbation Theory 

Quadratic perturbation Q 2 gives a finite shift in the first order PT. Denoting the 
unperturbed mode by \k x >= exp(ik x x), we have 

5 2 u 2 = c 2 < k x \Q 2 \k x >= - < V 2 > k 2 y . (11) 
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Averaging over angles gives a fractional frequency shift 

1 



s 2 = --M 2 , (12) 



where M 2 =< V 2 > /c 2 . 



The linear perturbation Q\ gives a non-zero frequency shift in the second order PT, 

r 2 2 V 1 ^ k x \Qi\k' >< k'\Qi\k x > . . 

^ =C L 12— 7^ • (13) 

For a sinusoidal shear flow with a wavenumber g, the sum is over k' — k x ± q, and the above 
expression gives 

^ - 4 < V- > ,|)^_±. (14) 



Sl= M 2 / o d n( l-a 2 ) ^ ^ , (15) 



Averaging this equation over the angle between k and q yields 

c 2\„2 u2 

where a is the cosine of the angle. For k > q/2, there is a pole due to the resonant scattering. 
The frequency shift is given by the principal value of the integral. The imaginary part gives 
the damping, which in the solar context corresponds to the line-width (Goldreich & Murray 
1994). We do not consider line- widths in this work. 

2.3. Results 

The integral is trivial in the long-wavelength limit, k q, 

si = (8/15)M 2 . (16) 

Summing the two contributions, Eq.(12) and Eq.(16), we get a total fractional frequency 
shift 

. = f Q M>. (17) 
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In the short- wavelength limit, k 3> q, Eq.(15) gives 

1 



si = ~M 2 : (18) 
and the total short- wavelength shift is given by 

s = ~\m 2 . (19) 

The fractional shifts in the asymptotic regions are wavevector independent. This 
is equivalent to the sound speed renormalization. The long-wavelength sound speed is 
higher than c, c 2 = c 2 + j| < V 2 >, and the short- wavelength sound speed is lower than 
c, c 2 = c 2 — | < V 2 >. These results were confirmed analytically (Appendix A) and by 
numerical calculations (Appendix B). 



3. THE SOLAR FREQUENCY SHIFTS 

The above results are not directly applicable to the problem of the solar p and /-mode 
frequency shifts due to interactions with turbulent convection. The Sun is stratified, and the 
solar convection is anisotropic and inhomogeneous. However the method of §2, second-order 
PT analysis of the density equation, is still applicable. 

In §3.1 we derive the basic density equation analogous to Eq.(8). In §3.2 we specify the 
model of the solar stratification in the upper convection zone. We give a general formula 
for the frequency shifts in §3.3, and calculate the f-mode shifts for a particular model of the 
solar convective turbulence in §3.4. 



3.1. The Basic Equation 



The aim is to calculate the eigenfrequencies of linear sound waves in a medium with 
an unperturbed density p , sound speed c, gravitational acceleration g in the presence of 
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a flow V. We use a planar model for which p and c depend on the depth x, whereas g 
is independent of the depth. The convective flow V satisfies the mass continuity equation 
V • (V Po ) = 0. 

Velocity perturbations v a and a density perturbation p (unlike §2.1 p is not normalized 
to p ) associated with the acoustic mode satisfy the linearized equations 

- iuv a + Vpdpv a + vpdpV a = — -d a c 2 p + —g a p, (20) 

Po Po 

- iup + dpVpp + d a p v a = 0. (21) 

Here and in what follows we use the convention that an operator acts on everything to the 
right of it, unless there is a bracket. We change variables, pov a — > v a , and define three 
operators 

F a = d a c 2 - g a , (22) 
G = d a V a , (23) 
H aP = S afi G + {d fi V a ). (24) 
Then the system (20), (21), can be written as 

- iuv a + H af3 vp + F a p = 0, (25) 

— iuop + Gp + d a v a = 0. (26) 
Now solve Eq.(25) for v in terms of p expanding up to the second order in V, 

v a = —F a p + H a/3 Fpp + HapH^Fyp. (27) 
iuj [lujy (icu) 6 

Substituting this expansion into Eq.(26), we get the basic perturbed density equation 
(analogous to Eq.(8)) 

QoP + Qip + Q2p = tu 2 p. (28) 
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Here the unperturbed operator is 

Qo = ~d a F a , (29) 



the linear perturbation operator is 



Qi = -d a H aP Fp - iuoG, (30) 
and the quadratic perturbation operator is 

Q 2 = —^daHapHp^F^. (31) 

Before we can use the PT technique on Eq. (28), we have to specify the stratification, 
find the unperturbed eigenmodes, and "Hermitize" Q . 



3.2. The Stratification in the Solar Convection Zone 

We use a polytropic model of the solar convection zone (see Goldreich & Kumar 1990). 
In this model the unperturbed (bare) sound speed is c 2 = gx/m, where x is the depth and 
m ~ 4 is the polytropic index in the hydrogen ionization zone. Then the unperturbed 
density equation, Eq. (28) with Qi = Q 2 = 0, reads 

- - (p" - k 2 p) + p' = —p, 32 

m m g 

for a horizontal wavenumber k. An eigenmode with n nodes in the vertical direction is of 
the form p^ = exp(— kx)x m ~ 1 L n (kx), where L n is the Laguerre polynomial of degree n. 
The dispersion law is u 2 , = (1 + 2n/m)gk. The modes with n > are the p-modes, and the 
n = mode is the fundamental mode (f-mode). 

The fundamental mode frequency does not depend on the value of the polytropic 
index m, oj 2 = gk. This a familiar dispersion law of the deep-water gravity waves. A 
remarkable fact is that the f-mode frequency is \fgk for an arbitrary stratification. The 
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solar stratification in the upper convection zone is not well known (Christensen-Dalsgaard 
1996), making interpretation of the p-mode shifts difficult. But the f-mode shifts must be 
mostly turbulent; we will concentrate on the f-mode shifts in §3.4. 

For a PT to be applicable, the unperturbed operator, Q , must be Hermitian. This is 
achieved by introducing the scalar product 

<Pi\P2>= J (frx l ' m p lP2 . (33) 

Knowing the unperturbed solutions and the scalar product we can now calculate the 
frequency shifts. 



3.3. Frequency Shifts 

We can use the unperturbed frequency in the perturbation operators Qi and Q 2 . The 
quadratic perturbation Q 2 gives a non-zero frequency shift in the first order PT. For a 
horizontal wavenumber k and a mode order n, we have 

5 2 uj 2 n {k) =< k, n\Q 2 \k, n > . (34) 

A non-zero frequency shift from Qi appears in the second order PT, 

j? 2/i \ <k,n\Q 1 \k',n' ><k',n'\Q 1 \k,n> . . 

wt) = £ ^m^m ■ (35) 

The total frequency shift is a sum of Eqs. (34) and (35). It remains to specify the flow V 
and do the integrals. 



3.4. Turbulent Shifts of the Solar Fundamental Mode 



We use a two dimensional model with a y-axis in the horizontal direction. For a single 
Fourier harmonic exp(iqy) of the turbulent flow V, the incompressibility of the mass flow, 
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d a (V a p ) = 0, is satisfied by the velocity field of the form 



V x = -iqV, 



(36) 



(37) 



x 



Since the sound speed scales as x 1//2 , we choose the stream function that scales as x 3 ^ 2 at 
small depths x, so as to make the Mach number of the turbulence finite. We suppose that 
for a given horizontal scale of the convection, only one row of convective cells exists in the 
vicinity of the surface. This means that \I/ has a maximum at the depth of the order of 1/q, 
and then goes to zero at larger depths. We choose 



where the y-dependence is omitted, [i measures the aspect ratio of the convective cells. We 
will sum the frequency shifts due to different Fourier components of the turbulence at the 
end of the calculation. 

For a polytropic model of the solar convection zone, the unperturbed eigenmodes 
are sums of products of exponentials and powers (PEP). We have chosen, Eq.(38), a 
convection stream function which is also a PEP. The perturbation theory requires a number 
of differentiations, multiplications and additions (which leave a PEP a PEP) , followed by 
an integration (a trivial algebraic operation for a PEP). We wrote a Numerical Algebra 
program which performs all these operations. The program was tested as follows. In the 
interval 0.1 < k/q < 0.6 the calculated frequency shifts coincided with the shifts calculated 
by the direct numerical simulations program described in Appendix B. Another reason to 
believe that we did not make a mistake, is that we succeeded in reconstructing analytical 
formulae with a correct algebraic structure from our Numerical Algebra program (Appendix 



We now present the results. For a single Fourier harmonic of the flow, the f-mode 



m = Ax zl2 e~ mx , 



(38) 



C). 
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fractional frequency shift is shown in Fig.l. The long- wavelength (k — > 0) shift is positive. 
The poles in the vicinity of k/q — 0.7 are due to resonant scattering. The resonant 
condition, u^(k ± q) — ^(k), is satisfied by two families of solutions. The first family is 

which in our case (m = 4) gives resonances at k/q = (n + 2)/(n + 4) = 0.5, 0.6, 0.67, 0.71, 
0.75, 0.78, 0.8, ... The first two of these are actually absent because the corresponding 
matrix elements vanish. The next five are seen in Fig.l, and higher resonances are not 
present because we summed the second-order PT shift, Eq.(35), only up to n! = 6. 

The second family of resonances is given by 



which in our case is k/q = 1 + 2/n = 3, 2, 1.5, ... These are shown in Fig. 2, where the 
"red" resonances are too weak to be seen. Fig. 2 also shows that the resonances disappear 
without a trace when integrated over a spectrum of turbulent modes. Even a rather narrow 
spectrum used in Fig. 2 reduces a complicated pattern of resonances to just two asymptotics 
and a transition region between them. 

The short-wavelength (high-degree) frequency shift does not depend on the aspect 
ratio of the turbulent eddies, /z, nor on the turbulence wavenumber, q. With an accuracy 
better than 1% , the high-degree shift is given by 



where M is the surface Mach number of the turbulence. The low-degree shift depends on jj, 
and q. It is given by the following exact formula 




(40) 



s 



uv 



— M 2 , 



(41) 



5 



M\k/q)\ 



(42) 



242/i 6 



-13- 



which is applicable for k/q <C /i. Smaller aspect ratios /i give higher maximal positive shifts 
s max , but we find that even the very small fj, — 0.1, gives s max ~ 0.1%, for a typical value of 
the Mach number, M — 0.1. This is smaller than the current observational accuracy. 

The smallness of the maximal positive shift may be explained as follows. At very small 
k, the mode's energy is concentrated at large depths, where the turbulence is too weak to 
produce a noticeable frequency shift. At larger k, the mode-turbulence interaction becomes 
stronger, but the negative short-wavelength shift takes over. 

The f-mode frequency shift is practically independent of the polytropic index m. The 
high-degree shift for m = 2 is s uv = — 1.1M 2 , which is close to the m = 4 result s uv = —M 2 . 
We also calculated the frequency shifts for the pi-mode and obtained similar results. 

The shape of the spectrum-averaged curve in Fig.2. seems to be correct. Shown in 
Fig. 3 are the data of Libbrecht, Woodard & Kaufmann (1990) and the theoretical curve of 
Fig.2 scaled to match the observations. In physical terms, the required scaling corresponds 
to the Mach number of the turbulent convection M ~ 0.3, and the turbulent length scale 
A ~ 3000A;m. 

4. CONCLUSION 

We have calculated the f-mode frequency shifts due to scattering of waves by turbulence 
using perturbation theory. A two dimensional model of convection was employed, but our 
methods are applicable to a three dimensional problem with only minor changes. Given 
a turbulent velocity field, one can use the perturbation theory to calculate the turbulent 
f-mode frequency shifts and line-widths. The main qualitative results for our simplified 
model are that the long-wavelength (wavelength larger than the turbulent eddy size) 
frequency shifts are positive but unobservably small and the short-wavelength fractional 
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shifts are negative and of the order of M 2 . 
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To test the homogeneous flow results, Eqs. (11), (14), consider sound wave propagation 
in a simple two-dimensional shear flow. Let the flow velocity along the y-axis be +V in 
the intervals 2n < x < 2n + 1, and —V in the intervals 2n — 1 < x < 2n, where n is an 
integer. Consider a sound wave with a wavenumber k propagating in the y direction. An 
exact dispersion relation can be obtained by usual methods. It reads 



A. A Simple Shear Flow 



ki tanh k\a k 2 tan k 2 a 



(Al) 



(uj-kV) 2 (uj + kV) 2 ' 
where a = 1/2 is the half- width of the layers, and k± t 2 are defined by 



{uo±kVf = c 2 (k 2 ±kl 1 ). 



(A2) 



These equations can be solved analytically in the second order in V, the result is 



3 

On the other hand, Eqs. (11), (14) predict 



5u 2 = V 2 (3k 2 - h 4 ). 



(A3) 



Since the average inverse squared wavenumber of the flow is 



5uj 2 = V 2 (3k 2 -4< 4> k A ). 



(A4) 



<-^>=-E(2^+l)^ 4 /E(2n + l) 



-2 



1 

12' 



(A5) 



the two results agree. 
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B. Numerical Simulations 

In two dimensions, for a single horizontal Fourier harmonic of the turbulent flow V, 
direct numerical solution of the eigensystem describing linear sound waves is possible. If we 
wish to determine the quadratic in M frequency shift of an acoustic wave with a horizontal 
wavenumber k, only two satellites of the wave, k — q and k + q (q is the turbulence 
wavenumber) must be included into the calculation. A system of ordinary differential 
equations results as described below. 

Write the system (25), (26) as 

P $ k = P q ^k- q + P- q ® k+q - (Bl) 

Here = (v x ,v y ,p) is the acoustic wave, the subscript k refers to the horizontal 
wavenumber. The unperturbed operator P does not depend on the flow V , Eq. (Bl) with a 
zero right hand side (RHS) gives an eigenmode in the absence of the flow. The perturbation 
operator P q is linear in V , the subscript q denotes the horizontal wavenumber of the flow 
V. To within the M 2 accuracy, the satellites satisfy 

P <S> k+q = P q $ k , (B2) 

P <£ k- q = P- q $k- (B3) 

Equations (B1)-(B3) were solved in three steps: (i) the RHS in (Bl) is zero, shoot to 
determine the unperturbed mode (ii) use the unperturbed mode $^ in the RHS of (B2), 
(B3), shoot to determine the satellites <&k±q- (hi) use both satellites in the RHS of (Bl), 
shoot to determine the corrected mode Solution obtained in (iii) should agree with the 
boundary conditions at the surface, x — 0, and at the bottom, x = H. 

In the absence of stratification (g = 0, c = const, p = const), for a number of various 
flows V, our numerical results agreed with the theory of §2. For the solar stratification, 
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reliable results were obtained only for 0.1 < k/q < 0.6. The reasons for a poor performance 
of the program outside of this interval seem to be as follows. Small k require too large 
depths H, because the mode is cut off at about A/k. Larger wavenumbers, k ~ q, are in the 
resonance region, and the resonances are sensitive to the position of our unphysical bottom 
H. However, numerical frequency shifts in the interval 0.1 < k/q < 0.6 were if-independent 
to a few percent accuracy. These results agreed with the exact (up to the numerical algebra) 
calculation used in §3. 



C. Exact Results 

Knowing the algebraic structure of the frequency shifts formulae, we used the Numerical 
Algebra program to obtain exact results. The simplest example is the first order PT shift 
5 2 c^o(^) gi ven by Eq.(34). This shift is an integral of a sum of products of exponentials and 
powers (PEP). Rescale the length, so that the turbulence wavenumber is q — 1, also take 
a unit-amplitude convection, A=l. Then the exponential function of the PEP under the 
integral (34) is e ~ 2 ( k +^) x . The highest power of x in the PEP is 7. Therefore, the resulting 
shift must be a polynomial in k and /i divided by (k + /i) 8 . The Numerical Algebra confirms 
this result. The reconstructed exact first order PT shift (34) is given by 

S2 = = 77 — , (a + axk + a 2 k 2 + a 3 k 3 + a 4 k 4 ), (CI) 
uj 4 (/i + ky 



where 



a = 20/i 2 , (C2) 

ai = -(20/i + 82/i 3 ), (C3) 

a 2 = -(250 + 191/i 2 ), (C4) 

a 3 = 44/i, (C5) 

a 4 = 363. (C6) 
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We used our Numerical Algebra Program to reconstruct analogous formulae for the second 
order PT sum, Eq.(35). These exact formulae were used to integrate out the resonances, as 
shown in Fig. 2. The matrix elements < n,k — q\Q i\0, k > are identically zero for k < q 
and n — 0, 1. This explains the absence of "red" resonances at k/q — 0.5 and at k/q = 0.6 
as seen in Fig.l. 

The lowest order in k long- wavelength shift is entirely due to the first order PT shift, 
Eq.(Cl). As it is, Eq. (CI) gives an asymptotic fractional shift s ir = 5k 4 //i 6 . The stream 
function amplitude A = 1 corresponds to the surface Mach number M = y / m(m + 3/2) = 11, 
for both of the Fourier harmonics q — ±1. This gives the mean squared Mach number, 
M 2 = 242, explaining the origin of 242 in the denominator of Eq. (42). 
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Fig. 1. — Fractional frequency shift (s) normalized by the turbulence Mach number squared 
(M 2 ) as a function of the wavenumber (k) normalized by the turbulence wavenumber (q). 

Fig. 2. — Resonances of the frequency shift are removed by averaging over a spectrum of 
turbulence wavenumbers. The spectrum was A 2 q ~ g 4 exp(— 2q 2 ) with a fixed aspect ratio 
H=l. 

Fig. 3. — Fractional frequency shift. The theoretical curve of Fig.2 was scaled to match the 
data of Libbrecht et al (1990) as explained in the text. 



